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Abstract. We introduce a quantum generalization of classical kinetic Ising 
models, described by a certain class of quantum many body master equations. 
Similarly to kinetic Ising models with detailed balance that are equivalent to 
certain Hamiltonian systems, our models reduce to a set of Hamiltonian systems 
determining the dynamics of the elements of the many body density matrix. The 
ground states of these Hamiltonians are well described by matrix product, or pair 
entangled projected states. We discuss critical properties of such Hamiltonians, 
as well as entanglement properties of their low energy states. 
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1. Introduction 

The paper that we have prepared for the special issue of the New Journal of Physics 
on " Quantum Information and Many-Body Theory" has two apparently independent 
motivations. First, it is motivated by the recent interest in many body quantum 
master equations and design of open systems for quantum state engineering and 
quantum simulations. Second, it concerns entanglement properties of quantum many 
body states corresponding to classical kinetic models with detailed balance and their 
generalizations to quantum master equations. 

1.1. Quantum master equations 

The quantum master equation (QME) is the basic theoretical tool to describe the 
evolution of open systems that undergo Markovian dynamics [1] [2] . While the initial 
studies of QMEs dealt with many body applications, over the last decades most 
research has focused on systems of few degrees of freedom, mainly because of the 
unquestionable complexity of the many body QMEs. In the recent fifteen years, 
however, there have been two waves of renewed interest in many body QMEs, related 
to the unprecedented progress of experimental control and engineering of ultracold 
atomic and molecular systems. On one hand, there has been a wealth of interest in 
the studies of various kind of cooling processes using QMEs. Zoller, Gardiner, and 
collaborators [3] studied in a series of papers the growth of a Bose- Einstein condensate 
in a trapped bosonic gas. One of us (M.L.), together with I. Cirac, P. Zoller, Y. Castin, 
and others, developed the theory of laser cooling of Bose gases [4], in particular in the 
so called festina lente limit [5]. Similar ideas were applied to the processes of laser 
cooling of Fermi gases [6] and sympathetic cooling [7J [8j |9] . 

On the other hand, several authors proposed to make use of the capabilities 
of modern quantum optics and atomic physics experimental methods to design and 
realize experimentally QMEs that allow to engineer interesting quantum many body 
pure states [10l HTJ |T2] . These pure states range from simple Bose Einstein condensates 
to the stabilizer states, MPS-PEPS states, or states with topological order. The first 
experiment realizing these ideas was conducted in the group of G. Rempe [13], who 
has been able to prepare a ID bosonic Tonks-Girardeau Q3] gas employing three- 
body losses. So far most of these new proposals concern ultracold atoms, molecules, 
in particular polar molecules, or ions, but recently they start to involve Rydberg 
atoms, that are particularly suitable for design and realization of 3-body, 4-body etc. 
interactions thanks to the Rydberg blockade mechanism. Weimer et al. [15] , proposed 
to use many-body quantum gates stroboscopically, employing the Rydberg blockade 
effect to engineer the topologically ordered ground state of the famous Kitaev toric 
code, the color code [16], and even to realize a quantum simulator of the U(l) lattice 
gauge theory. 

Motivated by these developments we propose in this paper a new class of 
many body QMEs generalizing classical kinetic models (in particular kinetic Ising 
models (KIMs); for a review see [T7]). Our QMEs have as stationary states thermal 
Boltzmann-Gibbs states of the underlying classical model. Note that these thermal 
states might correspond to very complex quantum states, for instance if the underlying 
classical model concerns kinetics of commuting stabilizer operators for the cluster 
states [IB] , as we discuss in Sec. [5] The diagonal elements of the density matrix in our 
models undergo dynamics equivalent to that of the underlying classical model. The 
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off-diagonal elements of the density matrix exhibit complex evolution, but fortunately 
can be grouped into independently evolving blocks. Similarly to kinetic Ising models 
with detailed balance, which are equivalent to certain Hamiltonian systems, our models 
reduce to a set of Hamiltonian systems determining the dynamics of the elements of the 
many body density matrix. In ID we identify classes of these models that are exactly 
soluble. In general, the ground states of these Hamiltonians are well described by the 
matrix product states (MPS), or pair entangled projected states (PEPS) [T9 ] l2"0" ] . 

It should be noticed that there exist previous efforts in the above direction. For 
instance, in Ref. [21] (see also Ref. [T7]) a general master equation for kinetic Ising 
model with an environment was derived. Under some assumptions, the diagonal 
elements of these formulations also reproduce Glauber's kinetic model, although no 
full solution of the equation was attempted. 

1.2. Entanglement in many body systems 

The studies of the role of entanglement in many body systems were definitely initiated 
by the seminal Ref. |22] , but they go back to the early works on area laws in quantum 
systems [23J (for a review see [Ml HH1 US] • Let us remind the readers that grounds 
states of (non-critical) many body systems with local Hamiltonians exhibit the area 
law. This means that if we divide the whole system into a subsystem A and the rest 
B where the size \A\ of A is large, but much smaller that that of \B\, and calculate the 
von Neumann entropy of the reduced density matrix of A, the latter will scale as the 
size of the boundary ("area") of A, Sa oc \dA\. The area law expresses the fact that 
away from criticality, correlations -and in particular quantum correlations responsible 
for entanglement- decay on short length scales. 

The situation is different at criticality, although so far only ID systems are 
fully understood. In one dimension, the area law at criticality may get logarithmic 
corrections so that the entropy of the block of size L scales as Sl oc clog(L), where 
the constant c can be related to the charge of the conformal field theory describing 
the corresponding critical behavior. In higher dimensions the situation is much less 
clear, and no universal laws at criticality are known [20]. Very recently, Masanes [26] 
proved under quite general conditions that for the ground states of systems with local 
Hamiltonians, the entropy of a block A is bounded from above by \A\ log A. 

There is a class of quantum states that always fulfill the area law in any dimension, 
despite the fact that they exhibit criticality for certain values of parameters [20] . These 
states are related to thermal states of classical Hamiltonians, such as Ising, or Potts 
models. In fact, any set of local quantum mechanical commuting operators can be 
used to build such models. For example, one can take stabilizer operators for cluster 
states [18] , or star and plaquette operators of the Kitaev model [16] . In the following 
we will mainly focus on Ising models, leaving the discussion of more complicated cases 
to further publications. 
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For Ising models, the states in question have the forrrjj] 

|*} = ^£exp[-/3tf(a)/2]|a), (1) 

where a — (a%, . . . , ojv) denotes a configuration of N Ising spins (|cr) stands for a 
vector representing a in the corresponding Hilbert space), H(o~) is the corresponding 
classical Hamiltonian, and Zn = Trexp[— (3H(o~)] is the partition function. As pointed 
out in Refs. [THEQ], and references therein, these states have the following properties: 

• They are associated to a family of classical kinetic Ising models (KIMs) that 
describe the approach to the classical thermal equilibrium state exp[— (3H(a)]/ Z^ 
and obey detailed balance. 

• The KIMs in question can be transformed by an appropriate ansatz into an 
equivalent problem of Hamiltonian dynamics in imaginary time (i.e. describing 
decay in time), with a certain quantum Hamiltonian H parametrically depending 
on (3. 

• The Hamiltonian H is non-negative and has one eigenvector given by the 
expression ([T]) corresponding to the eigenvalue zero. Away from criticality H 
is gapped, i.e. excited states have strictly positive energies. As f3 — > f3 c the 
gap vanishes, and at /3 C the Hamiltonian is gapless. The way the gap vanishes 
determines the dynamical critical exponent z of the associated KIM. 

• The ground state ([1} fulfills the area law and does not exhibit any special behavior 
at criticality. It can be exactly represented as a MPS in ID, or as a PEPS in 
higher dimensions [19] . 

Intrigued by the fact that the states ([TJ always fulfill the area law, we have 
attempted to look more closely into their properties, and the properties of their 
parent Hamiltonians H . One possible way is to study the entanglement properties 
of the excited states of H, but we have chosen another approach. We have generalized 
the KIMs to quantum models by defining a new class of QMEs, as explained in the 
previous subsection. These QMEs define new classes of parent Hamiltonians, and the 
grounds states of these Hamiltonians are expected to be again well described by matrix 
product or pair entangled projected states. This paper is devoted to the discussion 
of the critical properties of such Hamiltonians, as well as entanglement properties of 
their low energy states. 



1.3. Plan of the paper 

The paper is organized as follows. In Secion 2 we remind the readers of basics of KIMs 
with detailed balance, associated Hamiltonians, and more. In Section 3 we present 
the main result of this paper: we describe the generalization of KIMs to QMEs, and 

\ We would like to warn our readers against a possible misunderstanding: At issue are irreversible 
processes whose quantum behavior in general involves mixed states described by density operators. 
The generators in the pertinent master equations are non-Hermitian operators, precisely due to the 
irreversible dynamics. However, given detailed balance a master equation can be rewritten such that 
the generator becomes Hermitian and is then called a "Hamiltonian". Following widespread practice 
we shall indulge in calling a so transformed master equation a "Schrodinger equation in imaginary 
time" and talking about eigenstates of the " Hamiltonian" ; the resulting formalism then even allows 
to deal with superposition states, Eq. JTJ being the first example. Nevertheless, this paper is about 
density operators and damped motion. The "state" in Eq. |(TJ represents the canonical density 
operator. 
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we discuss the properties of the QMEs. In particular, we show how the equations for 
2 N x 2 N matrix elements of the density matrix split into 2 N equations of the KIM- 
type for functions of 2 N spin configurations. Detailed balance allows to transform 
these quasi-KIM equations to a Hamiltonian form, and we derive here the new classes 
of Hamiltonians. Properties of these Hamiltonians and their low energy states are 
discussed in Section 4. Here we show that Hamiltonians describing the evolution of the 
off-diagonal density matrix elements (coherences) are strictly positivc[§], which implies 
that the corresponding coherences decay to zero. At criticality we observe effects of 
critical slowing down of coherences; the Hamiltonians become gapless, implying that 
at criticality the decay is of the form of an exponential decay times algebraic tails. 
In some cases, criticality may even lead to survival of coherences for infinite times. 
In Section 4 we also discuss entanglement properties of the ground and excited states 
of the corresponding Hamiltonians. Our conclusions and outlook are contained in 
Section 5. 

2. Basics of kinetic Ising models 

As we discuss in Section 4, the methods developed in this paper can be 
straightforwardly generalized and applied to models other than Ising models: Potts 
models, classical clock models, or models employing commuting stabilizer operators. 
For simplicity and concreteness we will limit the discussion to Ising models, with Ising 
variables o~i = ±1, described quantum mechanically by the commuting Pauli matrices 
cr?. We will consider systems with a classical Hamiltonian H(a) following Markovian 
dynamics toward the thermal equilibrium state. The dynamics for a Markovian 
stochastic process is most conveniently formulated for the conditional probability 
P(cr, t\ao, 0) for the configuration a at time t, provided the initial configuration was 
(Jo at t = 0. This conditional probability allows to calculate all multi-time correlation 
functions of the process. Note that P(a, 0|<ro, 0) = 8 aaa: i.e. at the initial time 
conditional probability is obviously given by the Kronecker delta of the configurations 
a and o~o- In the following, in order to avoid too many arguments we will consider 
dynamics of probability of configurations P(a,t), which fulfills the same equation as 
P(o~, t\o~o, 0), but with a more general initial condition. 

2.1. Kinetic Ising models with detailed balance 

Here, we define in a more detailed way the classical kinetic Ising models and recall some 
of the literature results concerning particular examples of such models. In general, 
kinetic Ising models are defined by specifying the so called kinetic (master) equation 
for probability P(cr, t) of the form 

P(a, t) = £ M°' <r)iV> *) ~ °')P(^ *)] . ( 2 ) 

where the sum runs over all possible configurations a. The function w(a~ —¥ a'), 
hereafter also called transition probability, stands for the probability per unit time 
that the system changes its configuration from a to a 1 . 

§ An operator O S 3B{3%'') for some Hilbert space jriP is said to be positive if {tp\0\tp) > for any 
\il>) (z Then we say that O is strictly positive if (ip\0\ip) > for any |i/>) G J(f. For a Hermitian 
operator O acting on finite dimensional Hilbert space, both conditions are equivalent to nonnegativity 
and strict positivity of eigenvalues of O. 
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Usually one imposes the so called detailed balance condition (DBC) that says that 
at equilibrium the probability per unit time of a transition from a to a' is the same 
as the probability of transition in the opposite direction, a' — > a . Mathematically, it 
reads 

w(a' -> a)P eq (a') = w(a -> <i')P cq {a), (3) 

where P e q(a) is the equilibrium probability distribution P eq (a) = P(a,t — > oo). 

To support the above general remarks with some more detailed investigations 
let us now discuss some previously studied examples of such kinetic models. At 
the beginning let us focus on the one-dimensional Ising spin system, i.e., a chain 
of N spins uniformly distributed on a line. In this case a denotes one of the 2^ 
possible configurations of N spins and can be represented as a iV-dimensional vector 
(ax, ■ ■ ■ , <tn) with <7j = ±1 (i = 1, . . . , N). 

We restrict our attention to the case in which the behavior of ith spin is local, 
i.e., depends only on the nearest neighbors (generalizations to local models with 
next nearest neighbors are straightforward). We also assume that the probability 
distribution at equilibrium is 

P cq («) = ^-e-W') (4) 

with H denoting the classical (local) Ising Hamiltonian. In particular, we consider 
here the ferromagnetic Ising model which in ID corresponds to 

H(a) = (J>0). (5) 

i 

In this case the partition function has the explicit form Zn — 2 N (cosh N /3J + 
smh N p J) . The simplest possible process that may occur here is a single flip of the «th 
spin. Schematically this process can be stated as a — > Dia, where Di denotes the flip at 
ith position, DiO = (ax, • • • , — a,, . . . , ajsr). Also, let w(a — > -D,:cr) denote the transition 
probability for that process. The only processes that lead to the configuration a, 
appearing on the left-hand side of Eq. ([2]) are spin flips DiO — > a for any i = 1, . . . , N. 
The inverse type of processes can drag the system away from a. This means that the 
general kinetic equation ([2]) can be reduced in this case to a much simpler form 

N 

P(a, t)=^2 [w(Dia -> a)P(Dia, t) - w(a -> D l a)P(a 1 1)\ . (6) 
»=i 

The most general form of w(a — > Djcr) in the case that the interaction with both 
nearest neighbors is symmetric and leads the system to the equilibrium state (|4]) was 
shown to be [2"7] : 

w(a -t D t a) = T(l + 6ai-xa i+ x) [l - ^o^Oi-x + a i+l )] , (7) 

where 7 = tanh2/3J (notice that the value 7 = corresponds to infinite temperature, 
while 7 = 1 to zero temperature), \S\ < 1 and < T < 00. 

The case of S = was thoroughly investigated by Glauber, and shown to 
be solvable in the sense that all the relevant physical quantities can be computed 
analytically. In particular, the non-equilibrium expectation values and equilibrium 
correlation functions were determined. Moreover, this model was shown to have the 
dynamical exponent[j|| z = 2. Later on, the above model was treated in a series of 

I This means that the time dependent spin-spin correlation function decay on time scale as t^ ec 
behaves as t dcc oc £ z , where z is the dynamical critical exponent, and £ denotes the correlation 
length; § scales as (1 — 7)" - 1 / 2 when 7 — > 1. 
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papers [3 11 1521 125] for nonzero 6. In particular, it was shown in Refs. [23131] that the 
choice 6 — 7/(2 — 7) entails the dynamical exponent z^2. 

KIMs have also been studied in two and higher dimensions, although in these 
cases there is no known analytical solution. However, using efficient techniques 
precise numerical results have been obtained for relevant quantities such as the critical 
dynamical exponent |33) . 

2.2. Associated quantum Hamiltonians 

Interestingly, as it was show in e.g. Refs. [27J [Ml [3U] the detailed balance condition 
(J3J allows to rewrite the master equation ^ as a Schrodinger equation. In order 
to see that let us introduce the function 4>(a, t) related to the probability P(cr, t) 
through P(a, t) = ^/P cq (cr) <fr(a, t), where as above P eq is the equilibrium distribution. 
Inserting the latter in the master equation ([2]) and reorganizing slightly some terms, 
we arrive at 



(<J,t) = y [P^{(r)w{a' -► o-)pX(o-') 



,(</,*) 

= -yH„„,<t>(a',t) (8) 

a' 

with H denoting the real matrix with elements given by 

H aa , = y w(a" -> a 1 ) S aa , - P^ {a)w(a' -»■ a)P c \(o- r ). (9) 

a" 

Due to the detailed balance condition ([3]) one sees that H is symmetric and hence Eq. 
([8]) gives the aforementioned Schrodinger equation. 

Diagonalization of the Hamiltonian H is equivalent to the full solution of the 
corresponding master equation ([2]). Many of the previously discussed systems were 
investigated from this point of view. Below we recall some of the known results for 
the one-dimensional Glauber model. The Hamiltonian associated to Glauber's master 
equation ([5J with spin rates given by ([7]) has the form 

(7, S) = - Ty { [A( 7 , 5) - S( 7 , 5)aUal +1 ] of 

i 

-(\ + 8aUat +1 )[l - K«i - (10) 

where 7, 5, and T are specified as previously; a z and a x are the standard Pauli matrices 
and 

A{ 1 ,5) = — 1 , n -S, B 7,5) = 1-A(j,6). (11) 

For 6 = it was shown in Ref. [59] (see also Ref. [35 for another approach) that the 
standard procedure consisting of a Jordan- Wigner transformation followed by Fourier 
and Bogoliubov-Valatin transformations [371 13H] brings the Hamiltonian .ff (7, 0) to 
its diagonal form with non-interacting fermions. The eigenvalues are given by 

k 

M<lk) = yK, A„ = l-7Cos<z (12) 
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where q k denotes the ordered chain (qi < q 2 < ... < with each q t from 
{±ir/N,±3ir/N,...,±(N - l)n/N} for even N and {0, ±2w/N, ±4ir/N, . . . , tt} for 
odd N. The ground state, which is here the zero-energy state, is the one given by Eq. 
(Q~]). Moreover, the first nonzero eigenvalue is 1 — 7 and goes to zero for 7 — > 1 (zero 
temperature limit). 

The fermionic representation allows to go beyond Glauber's seminal work and 
solve exactly a whole class of master equations that are associated to the Hamiltonians 
given by Eq. (fT0|) . The inverse approach has also been explored: Given a quantum 
Hamiltonian that can be solved (at least partially), what are the corresponding 
classical master equations, and do they represent interesting physical systems |39j? 



3. Quantum kinetic Ising models 

Now we are in a position to proceed with the quantum generalization of the classical 
kinetic equations ©. Here we only discuss such a generalization for the Glauber 
master equation ([6]) with single-spin flips. However, further generalizations are 
obviously possible and are left for future work (see also Sec. [5]). We start by defining 
our notation. Naturally, as the computational basis in (C 2 )® N we take the eigenstates 
of TV- fold tensor product of a z . The basis consists of 2^ vectors hereafter denoted by 
\a) = I £7i, . . . , cttv) (c = 0, . . . , 2 N — 1) with er^ denoting an eigenvalue of of. After 
appropriate rescaling we may look at (ci, . . . , cat) &s the binary representation of the 
decimal number a. 

Let us consider the following master equation 

dg(t) ^ / 1 1 1 



dt 



£ {ofMOJ* 0(*)MO]* < ~ \ WO, <?(*)}} • (13) 



Here {•, •} denotes the anticommutator and Wi(cr z ) are quantum mechanical 
generalizations of spin rates introduced already in the preceding section^ (the Ising 
variables a = ±1 are replaced by the Pauli matrix a z ). Thus, all the spin rates Wi(a) 
(i = 1, . . . ,N) are diagonal in the computational basis \a). It is then clear that the 
diagonal part of Eq. (jTB"]) reproduce the kinetic equations ([U for all configurations a. 
Notice that Eq. (|T3| can be written as 

^=Y,{L ie {t)L\-\{L\L he m) (14) 

with the Lindblad operators given by Li = erf [wi(a z )} 1 ^ 2 . It is then clear that we 
deal only with the dissipative part of the general master equation describing Markov 
processes. 

Nevertheless, further generalizations to the full master equation are possible (see 
e.g. Refs. [21] and [17]) as, for instance, for the Ising Hamiltonian, the part governed 
by the Hamiltonian vanishes and such an equation would also reproduce the kinetic 
equations ©. 

In what follows we consider only the spin rates that satisfy detailed balance, i.e., 
w t {a z ) = itfi(D <£ r")exp[-2y3Jof (o?_! + £7? +1 )]. (15) 

% For the sakes of clarity and simplicity in the case of quantum spin rates we change sightly the 
notation from w(cr — ► Dicr) to u>;(ct z ). According to this convention ui(D;cr — > cr) is replaced with 
Wi(DiO z ). 
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Let us now proceed with solving the above master equation. As we will see below, 
our method relies - as in the case of classical kinetic equations - on the observation 
that the whole equation can be brought to a set of 2 N Schrodinger equations. For this 
purpose it is convenient to use the isomorphism between matrices from M^(C) and 
vectors from C d ® C d . More precisely, we can represent the density matrix g{t) as 

g(t) = J2Mt)U«W)(°\ \Q(t)) = ^>Wk ? H|a). (16) 

cr.er (T,a 

This form of g(t) allows to rewrite the master equation (|13|) as the following matrix 
equation 

lei*)) = E I ? 5 ? VMV)*»i<?*) ~ \ [Wi{a z ) + Wi(a z )]\ \g(t)), 

(17) 

It should be emphasized that operators corresponding to the "tilded" and "nontilded" 
spins act on the "tilded" and "nontilded" kets in the vectors representation of the 
density matrix git) P^|) (for instance afaf\s)\s) — (Xi\s)af\s)). 

It is evident that the matrix appearing on the right-hand side of Eq. (jTTJ) is not 
Hermitian. In order to bring it to Hermitian form we can use the detailed balance 
condition (TTSl) . This suggests the following transformation (see e.g. Ref. |29j ) 

\g(t)) = a&-(P/4)[H(a)+H(?)]}W)), (18) 
with % denoting the quantum generalization of the Ising Hamiltonian given by Eq. 
([5]). Application of the above to Eq. (fT7f leads us to 

i#)> = E U-v-[vr{° z )\H<y)] h - + m° z )}} \m), 

i ^ * 

(19) 

where Vi(a z ) — Wi(a z ) exp[f3Jcr z (a z _ 1 + of +1 )]. Now, due to the detailed balance 
condition, it is clear that af and Vi(a z ) commute, and consequently we have a 
Schrodinger equation \ip(t)) = —H\ip(t)) with Hermitian H. Let us notice that the 
above procedure, being just a "quantum" generalization of the procedure described 
in Sec. I2.2| replaces the problem of solving the master equation for N spins with the 
problem of solving the corresponding Schrodinger equation for 2N spins. 

As we will see below the above equation can be further simplified. Namely, 
identifying operators that commute with H we can split the equation to a group of 
2 N Schrodinger equations. For this aim let us notice that H commutes with ofo z 
(i = 1, . . . , N). We can therefore introduce new variables Tj = of of (i = 1, . . . , N) 
which are constants of the motion. Then, one sees that the tilded variables can be 
expressed by r = (ti, . . . , tjv) and a variables as af = TiO~\ for any i. In other words 
we have replaced a and a by r and cr, of which r is conserved. The Hamiltonian H T 
then takes the form 

Hr= -E{^^(^)] 4 K(^ Z )]*-^K-K)+^(^)]}: (20) 

where to~ z denotes TiO~ z (i = 1,...,N). Let us make a comment on the notation. 
One sees that there are 2 N different configurations of r-spins. We label them by the 
natural numbers 0, . . . , 2 N — 1; since we shall repeatedly have to refer to the two fully 
homogeneous configurations we simply denote these as r = (all r-spins up) and 
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t = 2 N — 1 (all T-spins down) , without further specifying the association of naturals 
to configurations. 

By definition the r = configuration represents the equal values of a and a 
spins. Consequently, via Eq. (1161) . the corresponding Schrodinger equation describes 
the diagonal elements of the master equation (fT3]) and thus recovers the classical 
kinetic equations ([6]). In the remaining cases of r configurations the er and er variables 
have to differ on at least one position implying that the Schrodinger equation related 
to any r ^ describes the set of 2 N off-diagonal elements of g(t). 

Thus, by identifying TV operators commuting with H we brought the solution of 
the general master equation (|13l) to the problem of diagonalization of 2 N Hamiltonians, 
each of dimension 2 N x 2 N . 

Let us now concentrate on a particular choice for the transition probabilities 
Wi(o~ z ). In what follows we investigate the quantum version of the rates given by ([7]) 
(oi are replaced with erf ), that is 

Wl (a z ) = T(l + 5a*_ l( r* +1 )[l - KW-i + ( 21 ) 
with r, 5, and 7 defined as before. Putting this to Eq. ([20"]) and after some algebra we 
get 

# T ( 7 , S) = - TJ2 { [Ml, S) - Bifr, 5)al 

i 

-Sfin-m+^al^} , (22) 



2 -(i + 5)o-l [Ht^o-U + f{nn +1 )at +1 ] 



where 



and 



MciiS), n-x = n+i, 

Ml,S)={ (23) 



B(*y,6), rt-i = r i+ i, 
LM-.rf) : { (24) 

0, n-i = -n +1 . 

with A(j,S) and 5(7,(5) defined as in Eq. (JTTJ) and /(a;) = (1/2) (1 + x). One 
sees that for r = or r = 2 N — 1 the function / equals one for i = 1, . . . , N and 
therefore we obtain the Hamiltonian given in Eq. (fpQ|) . Let us stress again that 
the Schrodinger equation corresponding to r = describes the diagonal elements of 
the master equation (IT31 and thus gives the classical kinetic equations ©■ All the 
remaining rs correspond to the off-diagonal elements of g(t). 

For the interesting special case 6 = all Hamiltonians H T can be diagonalized 
using the results of Ref. [JT]. This is because for any r = 0, . . . , 2 N — 1, the Jordan- 
Wigner transformation [3 6) brings H T (-y,5) to Hamiltonians which are quadratic in 
the fermion operators 



i-l 



n 



|< : at = l(aV + ial). (25) 



\ j=1 

This, by virtue of the results of [51], means that diagonalization of any of £f T (7,0) 
reduces to the diagonalization of an N x N matrix. Since, in principle, the latter can 
be performed numerically efficiently, the Hamiltonians (|22[) are diagonalizable. 
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For the case of <5 7^ the Jordan- Wigner transformation gives us Hamiltonians 
containing terms which are quartic in Cj. In this case we can try the transformation 
proposed by Siggia [35] ■ For this aim we introduce new bond variables Zj = ofof^, 
erf = Xi-iXi, and Yi = — L?jJQ. They commute when their indices are different 
(different bonds) whereas for coinciding indices (same bond) they obey the algebra 
of Pauli matrices (Xf = 1, XjYi = iZi etc). The bond variables allow to rewrite the 
Hamiltonian from Eq. (|22|) as 



Generally this is the anisotropic Heisenberg Hamiltonian with some external field. 
In some particular cases such Hamiltonians are analytically diagonalizable (for a 
pedagogical review on the Heisenberg model and the Bethe ansatz that is used to solve 
it see Ref. [H]). It seems that the general case 5^0 defies exact diagonalization. 
However, in the following subsection we show that for 6 = — 1, H T (j,S) can be 
diagonalized analytically for all values of t. 

4. Properties of the Hamiltonians associated to the QME 

Here we study some of the properties of the Hamiltonians H T . First, we show that all 
of them are positive and in a majority of cases even strictly positive. We also discuss 
the cases in which H T have zero eigenvalues and thus identify possible stationary 
states of the master equation (|13l) . Then we study fully soluble case of the so- 
called energy conserving spin flips, where all Hamiltonians H T t = 0, . . . ,2 N — 1 
are diagonalizable. Finally, we numerically investigate properties of entropy of ground 
states of our Hamiltonians. 

Jf.,1. Zero eigenvalues of H T - stationary states of the evolution 

It is interesting to ask which states survive the evolution, that is, which of the 
off-diagonal elements of g(t) do not decay to zero with t — > 00. This can be 
done by studying the positivity properties of H T and in particular their eigenstates 
corresponding to zero eigenvalues. Generically, as we will see below, for any r except 
for r = and t = 2 — 1, the Hamiltonians H T are strictly positive except for two 
different (7, <5)-points, namely, (5 = 7 = and 5 = 7 = 1. On the other hand, for 
r = or t = 2 N — 1 it is known 29, 28 that the corresponding Hamiltonian has zero 
eigenvalues (two-fold degenerate) for any value of 7 and i5. 

Let us treat the two cases r / and r 7^ 2 N — 1 in a more detailed way. For 
this purpose let us notice that we may write H T (j,6) as H T (^,S) — Hr ^ (7, S) 
and study positivity of each Ht(j, 8). On the other hand, for any r 7^ 0, 2 N — 1 we 
may divide all such terms into two groups, the one consisting of (7, 8) for which 
Ti-i — t,+i and the one for which Tj_i = — t»+i. More precisely, we can write H T (7, 8) 
as 




[l- 1 (l+5)f(r i - in )Z i }}. 



(26) 



^ =) (7,^) + 



]T H^(j,8). (27) 



{i\ T i-l= T i+l} 



{i|Ti_l = -T i+ l} 
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In the first case of equal rs we have 

H?>=\<y,S) = T{1 - 2(1 + «)/(r i _ir i )«7f + af +1 ) + Sa^ +1 

-[A i ( 1 ,6)-B i ( 1 ,S)at_ 1 a! +1 ]at}. (28) 

One finds that Hr\j, 5) has four different eigenvalues, each two-fold degenerate. For 
— 1 < 6 < 1 their explicit forms are Ai = 0, A2 = 2r(l — 6), 



xij = r(i + 5) 



(29) 



For the chosen parameter region A2 > 0, while A3 ' is manifestly positive. On the 
other hand, to prove nonnegativity of A4 it suffices to notice that 1 — 7 + [(1 — 
r i „ir l )/2] 2 7 2 < 1 - 7 + 7 2 < 1 for 7 < 1. 
In the case of r,-_i = — Ti+i we have 

J ff(^)( 7 ,<5) = r {1 - 5(1 + *)o? (/(r^T^t! + /(nn+xX+O 

-Vl-52 $/l^fa?}. (30) 

Here one finds that since t,;_i = — Tj+i, one of the factors /(rj-iTj) or f(nTi + i) 
vanishes. Without any loss of generality let us assume that f(nri+i) — 0. Then, 
obviously f(ri-\Ti) = 1 and the above Hamiltonian can be brought to 

H^( 7 ,S) = 1 - i(l + 6)a?_^l - i/T^af. (31) 

It has two different eigenvalues (each two-fold degenerate) of the form 



1 



A± = 1 ± - Y 7 2 (1 + S) 2 + 4vl^7 2 (l — S 2 ). (32) 

One sees that obviously A^' > 0. To see that also A_ is nonnegative it suffices to 
notice that for < 7 < 1 and \S\ < 1 the maximal value of the function of 7 and 
S appearing under the sign of square root is four and is this value is attained for 
7 = 8 = and 5 = 7 = 1. This also means that 7 = 5 = and 7 = 5=1 are the only 
points for which Hr (7, 6) can have zero eigenvalues. 

In conclusion, it follows from the above analysis that all (7, S)s are positive 
and therefore the Hamiltonian -ff r (7, 5) is positive for any r = 0, ...,2 — 1, and 
parameter region specified by the conditions — 1<<5<1,0<7<1, and T > 0. 
Moreover, it follows that for all r ^ 0,2^—1, the Hamiltonians i? r (7, <5) are in 
general strictly positive except the points (5 = 7 = and 6 = 7 = 1. Thus these 
two points together with two values of r = 0, 2 N — 1 are the only possible cases when 
H T (j, S) can have zero eigenvalues. Let us discuss shortly each of these cases. 

For r = and r = 2 N — 1, as previously noticed, the corresponding Hamiltonian 
(|10[) has a zero-eigenvalue eigenstate given by Eq. ([IJ. The case r = (t = 2 N — 1) 
corresponds to the diagonal (anti-diagonal) elements of g(t), which as it follows from 
Eqs. ([T]) and (|18|) are of the form (taking into account the normalization) 

\Qr(t)) = ±- e"^ (ff) £ - 0, 2 W - 1). (33) 

Z N ^ 

For the zero temperature limit (/? — » 00), the above becomes the well-known ./V-partite 
Schrodinger cat state (or Greenberger-Horne-Zeilinger state) 

|^) = i=(|t)^ + ||)^) (34) 
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with | f) and | i) denoting the eigenvectors of a z corresponding to the positive and 
negative eigenvalue, respectively. Moreover, for 7 = 1 the two-fold degeneracy appears 
in the ground state of Hi(l,S) and the second zero-energy eigenstate is 

i^) = ^(itr-nr)- (35) 

Since generically the remaining off-diagonal elements of g(t) vanish in the t — > 00 limit 
(except for the already mentioned values of 7 and 5), we have an example of a state 
that is a stationary state of the master equation and becomes genuine multipartite 
entangled state in the limit of zero temperature. 

Let us now treat the cases of 7 = 6 = and 5 = 7 = 1. In the first one, the 
Hamiltonians if T (7, S) simplify significantly and are of the form 

iT T (0,0)=r^(l-a?) (t = 0,...,2 w -1), (36) 

i 

meaning that the only zero energy state is | — >-)® iv where | — >) and | 4—) stand for 
the normalized eigenstates of a x corresponding to positive and negative eigenvalue, 
respectively. Since, as it follows from Eq. (fT8|) . for 7 = 1 it holds that \g T (t)) = \ipt{t)) 
for any r and therefore \g T (t)) = (l/VW)] ->)® w for r = 0, . . . , 2^-1. Consequently, 
the stationary state is fully separable and is given by g s t{t) = with -P|_>) denoting 
a projector onto |— >) . 

The case of S = 7 = 1 is a little bit more difficult. Now, from Eq. (|2"21 it follows 

that 

H T (1, 1) = r£ {1 - of [/(Ti-iTiXa + f{T i T i+1 )at +1 ] 

i 

+/(r J _ 1 r l+1 K_ 1 < +1 ]}. (37) 

It is clear that this Hamiltonian is diagonal in the standard basis |<r) in (C 2 )®^ and 
thus we can look for the eigenstates among the standard basis in (C 2 )®^. Interestingly, 
using the previously introduced bond variables Zi = o-f_ 1 crf, in the case of periodic 
boundary conditions this Hamiltonian can be brought to the antiferromagnetic Ising 
Hamiltonian with magnetic field 

H T (1, 1) = TN + fin-iTi+xjZi-tZi - 2Tj2f(T^in)Zi. (38) 

i i 

Let us concentrate on the zero-energy eigenstates of H T (1, 1). The latter can be 
found by solving the corresponding equation for eigenvalues of a z . This, however, due 
to the fact that in this equation each term under the sum is positive, means solving 
of the following set of equations 

f{Ti-\Ti)ai-\ai + f(TiTi + i)aia i+ i — /(Ti_iTj + i)<7i_i0i + i = 1 

(i = l,...,N), (39) 

where Oi stands for the eigenvalue of of . For instance, for all configurations of r that 
consist of blocks of length greater or equal two separated by the domain walls, one 
of the possible solutions is given by o~i — r^. This is because in such case the above 
set becomes f(n-iTi) + f(nT i+ i) - /(t,_it j+1 ) = 1 (i = 1, . . . , N). The only possible 
triples (T,_i,Ti,Tj+i) that can appear in the discussed case are ttti Hi, Hi, itt, 
Hi, and H^. It is clear that for all of them these equations are satisfied. One may 
also easily check that generically the zero-energy eigenstates are degenerated. As a 
result it seems that in the case of 7 = S = 1 there is a variety of states g(t) that are 
stationary states of our master equation. 
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4-2. The case of energy conserving spin flips 

Here we consider the case of energy-conserving spin flips, that is, when 5 = — 1. One 
sees now that the Hamiltonians H T simplify significantly and read 

ILrfr, -1) = TN - fi {Xi-xXi + Yi-xYi + , (40) 

i 

where f[ = /(rj-iTi+i). For r = and t = 2 N — 1 this is the isotropic ferromagnetic 
Heisenberg model with spectrum shifted by TN. In the general case, however, 
H T (p/,—l) depends on the numbers fi which are either zero or one depending on 
the configuration r. To deal with this it suffices to notice that for any configuration 
of t 7^ 0, 2 N — 1 some set of indices {ii, . . . , ik] exists for which fc k = and between 
these zeros the fi are constant and equal to one. It is clear from Eq. P0|) that such 
zeros divide the Hamiltonian into a sum of "smaller" commuting Hamiltonians. More 
precisely if = for j = 1, . . . , k then 

21 — 1 12 — 1 

H T ( 7 ,-l) = TN-T^2S n .i-S n -T S n -!'S n -... 

ri— 1 n=i\-\-l 

- J2 S n-1-Sn, (41) 
n=ife_i+l 

where S n — [X n , Y n , Z n \. It follows that S n -i-S n commutes with S m -i-S m whenever 
\n -m\>2. 

This means that for a given configuration of r we can split H T (j, —1) into a group 
of commuting isotropic ferromagnetic Heisenberg Hamiltonians with free ends. Such 
Hamiltonians can be treated using the so-called Bethe ansatz [33]. To visualize what 
we have just said let us consider the following illustrative example. Let us assume the 
periodic boundary conditions in Eq. (|4U|) and let the r configuration together with 
the corresponding chain f[ be given by 

(42) 

Then the Hamiltonian becomes (forgetting about the constant part and putting T = 1) 

il4 

H=13 ' 



T 


t t t 


I 


1 


t 


t t 1 


t I t t t 


fi 


1 1 











1 1 


110 11 



i?r(7: -1) = -QI i= i Si-1 ■ Si + S 6 ■ S 7 + Y^ltg S i-1 ' S i + Y,i=13 S i-1 ' S i)- 



4-3. Entropy of the ground state of the off-diagonal Hamiltonians 

Here we take a brief detour and study the entanglement between parts of the ground 
state of the Hamiltonians H T (j, S) that control the dynamics of the off-diagonal terms 
of the QME. In order to study our system for arbitrary values of S, we take advantage 
of the matrix product state structure that the ground state of these Hamiltonians 
has. For this, we extended the time evolving block decimation (TEBD) algorithm of 
Vidal 44 to include next-nearest neighbor interaction. We then performed a variable 
step imaginary time evolution to find the ground state of open boundary chains of 
length N large enough so that the results become independent of size. Our goal is 
to compute the bipartite entropy S = tr log 2 pl, where pl is the reduced density 
matrix obtained after tracing out N — L neighboring spins of the chain. For this, 
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instead of obtaining and diagonalizing p^, we make use of the Schmidt coefficients 
that appear explicitly in the TEBD representation of the statcQ 

S 
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l 


r- 




H 








iii — — 




« 
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Figure 1. Bipartite entropy S(L) = tr p^ log 2 Pl, where p^ is the reduced 
density matrix obtained after tracing out N — L contiguous spins of the chain, 
for the Hamiltonian d 2 2 1) in the 5 = case, and where r is a configuration 
where all but two components are equal. The curves, in ascending order, are 
for 7 = 0.1,0.3,0.5,0.7, and 0.9. The entropy reaches 1 for criticality at 7 = 1. 
The system is a chain of 90 spins. Notice how at the point of the flipped rj's the 
bipartite entropy is reduced but does not go to zero. In thin lines, for comparison, 
is the same entropy but for a configuration with all tj 's equal to one (the classical 
KIM, here the diagonal part of the QME). 

We study two representative cases of r configurations, one in which only two rj's 
are different than the rest, and one where half neighboring r^'s are equal between 
them and different than the other half - a case that we call a " domain wall" . We also 
concentrate on two relevant values of 5, first 5 = (the Glauber model), shown in Figs. 
(1) and (2), and second, the temperature dependent 6 = 7/(2 — 7) (the Haake-Thol 
model that gives the dynamical critical exponent 2 = 4), shown in Figs. (3) and (4). 
For comparison, we also compute the bipartite entropy of the diagonal component of 
the QME, i.e., the classical KIM. 

We observe that the maximum value of entropy grows with 7 and approaches 
unity for criticality (7 = 1). The block length L at which the entropy saturates also 
appears to be rather small, about 5 sites, although we expect this to grow near the 
critical point. 

In all cases of S a domain wall appears to de-entangle the two parts of the chain, 

+ Any pure state \ip) of a quantum system partitioned into two parts A and B can be written in 

its Schmidt decomposition form, \ip) = X/X=l A^ 2 \4'a^' 1 ) \<Pa ^)-> where {It/'!' 4 ')} and {|v'q S '}} arc 
two orthonormal bases of the Hilbert space of parts A and B respectively, \ a are non-negative real 
numbers (the Schmidt coefficients) , and \ i s the smallest of the dimensions of the Hilbert spaces of 
A and B. If part B is then traced out, the reduced density matrix of part A can be written as a 
diagonal matrix in the {J^! 4 ')} basis, with the Schmidt coefficients in the diagonal. Since the TEBD 
algorithm basically stores the Schmidt coefficients and the Schmidt bases for all possible bi-partitions 
of the system, then computing the von-Neumann entropy of the reduced density matrix p^ for any 
L-value is quite easy. 
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Figure 2. The same as Fig. 1 but for a configuration of Tj's where half are 1 
and the other half are — 1 . Notice how the entropy does go to zero at the domain 
wall. 
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20 40 60 80 



Figure 3. The same system as Fig. 1, but for 8 = 7/(2 — 7), where the classical 
KIM shows an anomalous dynamical exponent 2 = 4. In this case, the flipped 
T{ induces a spike in entropy near the domain wall, indicating some interesting 
quantum correlation existing between domain walls. Still the curves, in ascending 
order, are for 7 = 0.1,0.3,0.5,0.7, and 0.9. In thin lines, for comparison, is the 
same system but for a configuration with all r^'s equal to one (the classical KIM, 
here the diagonal part of the QME). 
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giving zero entropy for a partition right at the domain boundary. On the other hand, 
a single flipped shows some residual entropy. Therefore, the separation of the 
Hamiltonian into commuting parts shown above for 5 = — 1 is not possible for general 
t configurations. 

S 

0.8 r 




20 40 60 80 



Figure 4. The same system as Fig. 3, but for a domain wall configuration for 
r. In contrast to the case shown in Fig. 2, here the entropy grows sharply before 
the end of the domain. Nevertheless, it goes to zero exactly at the domain wall, 
again indicating that in this case the domains are not entangled. 

The case with an anomalous dynamical exponent z = 4, S = 7/(2 — 7), shows 
some interesting behavior of the entropy near the domain walls or the flipped Tj . In 
particular, we observe an entropy spike before the end of the domain. The effect is 
acutely pronounced near criticality. 

5. Conclusions and Outlook 

We have presented here a novel class of QMEs that have the following properties (cf. 
e.g. Refs. [21] and [17]): 

• The diagonal elements of the density matrix in the " computational" basis follow 
dynamics of a certain classical kinetic model. 

• The dynamics of the off-diagonal matrix elements splits into blocks, described by 
"kinetic-like" models. 

• For models fulfilling DBC, the dynamics can be transformed into a Hamiltonian 
dynamics in imaginary time. 

• The ground and low excited states of the resulting Hamiltonians fulfill area law, 
and can be well described by MPS or PEPS methods. 

Our results suggest several directions of investigations, which we would like to 
follow in the future: 

• In the present paper we have focused on Ising spins, and on generalized QME 
associated with KIMs. Generalizations to models associated with kinetics of more 
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general set of commuting operators, such as stabilizer operators, are possible and 
interesting. 

Several presented models admit exact solutions via Wigner- Jordan transformation 
or Bethe ansatz a la [351 129] . and/or approximate treatment using variational 
methods a la [25] . These methods should allow for more rigorous analysis of the 
novel types of Hamiltonians. 

Especially interesting are the two-spin-flip models such as those that conserve 
the magnetization [40 or energy |30j . In particular, for the energy conserving 
model in ID the generalized QME reads: 



V 71 - «+l 

(43) 



-- {l-^i-i^i<+i^i +2 ^(i)} 



The high degeneracy of the ground state in this model allows to expect the 
appearance of long-living coherences. 

• Generalization to models that do not fulfill DBC such as exclusion models (see 
Refs. [45]|46]) is possible. An example of somewhat "hidden" DBC is the QME 
of the form 



I— 1 



(44) 



where — l<a<lisa free parameter. This model corresponds to ID anisotropic 
Heisenberg model. 

Last, but not least, physical implementations of the considered models with 
ultracold atoms in optical lattices, or ions in trap arrays, or Rydberg atoms are 
feasible, and will be studied. 
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